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Abstract 

The stability properties of models of spontaneous mirror symmetry breaking in chemistry are 
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solutions. This parameter is identified in each case, and the constraints on the chemical rate 
constants determined from dynamic stability are derived. 
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I. INTRODUCTION 

Ever since the Soai reaction, first published some years ago [1], there has been an ever 
increasing number of researchers inspired by this important landmark experiment, who have 
been searching for other reactions capable of exhibiting the spontaneous emergence of chiral 
asymmetry in closed systems. At present, there is a gradual and mounting experimental 



evidence in favor of this 
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showing that absolute asymmetric synthesis [5|| is possible 
for autocatalytic reaction networks in closed mass reacting systems and even suggests that 
chiral amplification could occur in chemical equilibrium scenarios [6]. On the theoretical 
side, kinetic schemes which are extensions of the elemental Frank model [7|, are able to 
reproduce the main features of the mirror symmetry breaking behavior of the Soai reaction 



The key ingredients of theoretical models of mirror-symmetry breaking processes in chem- 
istry j^l include reactions in which the chiral products serve as catalysts to produce more 
of themselves while inhibiting the production of their enantiomer or mirror-image counter- 
parts. The two basic models we analyze here differ in the way this inhibition is produced. 
Frank's original model 0,0 , [h], [if], [jjj], involves the autocatalysis of the two enantiomers, 



denoted herewith as L and D, and mutual inhibition or antagonistic effects between the two 
chiral species. This mutual inhibition occurs through the formation of heterodimers that 



are removed from the reacting system. In the case of limited enantioselectivity [13], the 
required antagonism arises through the recycling of an enantiomeric pair of monomers back 
to a single chiral monomer and a prochiral substrate molecule. 

The purpose of this paper is to determine the stability of the steady states for vari- 
ous reaction schemes that have been proposed as models for mirror symmetry breaking in 
chemistry. This requires knowledge of only the static solutions of the kinetic equations. 
Then, differential equations for the arbitrary time-dependent fluctuations about these final 
static solutions are straightforward to derive to first order in the fluctuations. These time 
dependent fluctuation equations are expressed in matrix form, and the eigenvalues of the 
corresponding matrix, evaluated on the static solutions obtained previously, signal unam- 
biguously the stability (or instability) of the racemic or chiral solution under study. Here 
we are interested in understanding the chemical factors controlling the symmetry break- 
ing transition from a racemic (equal proportions of the left and right-handed molecules) to 
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chiral states, the latter characterized by unequal proportions of the two enantiomers. For 
each model considered, we identify the critical parameter that controls the transition from 
racemic to chiral states. These dimensionless parameters are expressed as simple ratios of 
the chemical reaction rates. When the control parameter falls below a critical value, then 
chiral amplification results. In the latter case, the racemic state is unstable, and an initial 
small perturbation is sufficient to tip the system over into one of its two equally likely sta- 
ble chiral states. Exactly equal proportions of the two chiral enantiomers never occur in 
practice, and this inevitable chiral imbalance yields an initial statistical enantiomeric excess 
14j . Moreoever, in imperfectly mixed spatially extended systems, diffusion-limited noise is 
intrinsic to the reacting and diffusing system itself and is sufficient to drive the symmetry 

nnn 

breaking [15|, [16|, 117(1 . In either case, the noise need not be put in a-posteriori. 

First the study of open systems at constant concentration of the achiral reactants is 
carried out. Open systems provide a clear insight to the chemical conditions necessary to 
achieve spontaneous mirror symmetry breaking; as we will show, the critical parameter for 
the bifurcation behavior in a closed system is concentration dependent. The calculation 
of the eigenvalues in the closed system case are generally more difficult to obtain and in 
many cases, they cannot be expressed in a manageable analytic form. The reaction schemes 
analyzed below correspond in the most part to reversible chemical reactions, in contrast to 
the irreversible Soai reaction. The results obtained here should be of practical value to aid 
in the search for new reactions exhibiting chiral symmetry breaking which might eventually 
provide valuable insight to the much more difficult problem of the origin of biological chirality 
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The importance of cyclic catalytic reactions as conditions for life is elegantly argued 



19]. 



This paper is organized as follows. After introducing the general reaction scheme in 
Sec [Til we consider open flow systems, and use the constancy of the prochiral substrate 
to cast the kinetic equations in dimensionless form in Sec IIHI This has the advantage of 
reducing the number of parameters and allows for all results to be expressed in terms of 
ratios of the rate constants. This is a most welcome feature when calculating solutions and 
the corresponding eigenvalues of the stability matrix. The stability properties of a sequence 
of Frank-type models followed by models with limited enantioselectivity are presented in 
Sec IIVI The symmetry breaking in all these models yields the classic bifurcation diagram 
which we illustrate for the limited enantioselectivity model. For closed systems, a different 
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rescaling of time and the concentrations is needed to obtain dimensionless kinetic equations. 
The constant mass constraint makes the algebraic analysis in this case much more involved, 
and a few particularizations of the general reaction scheme are treated in Sec El Conclusions 
are drawn in Sec IVII 



II. THE GENERAL KINETIC MODEL 



The general system we will study for spontaneous mirror symmetry breaking is that 
defined by Eqs. ([HIS]) below and it implies the following four pairs of reactions: a straight non- 
catalyzed reaction Eq.([T]), enantioselective autocatalysis Eq.(j2]), and a non-enantioselective 
autocatalysis Eq.([3]), where A is a prochiral starting product, and L and D are the two 
enantiomers of the chiral product. We also assume reversible homo- and heterodimerization 
steps in Eqs. ( MoT) ; L 2 ,D 2 represent the two chiral enantiomeric homodimers and LD is the 
diastereomeric, achiral heterodimer. The ki denote the reaction rate constants. 
Production of chiral compound: 

fel ki 

A Z L - A C D - (l) 



Autocatalytic production: 



&2 &2 



L + A "~ L + L, D + A " D + D. (2) 

k-2 k-2 



Limited enantioselectivity: 



^3 ^3 



L + A," L + D, D + A 7~ L + D. (3) 

K-3 K-3 



Homo-dimerizations: 



Hetero-dimerization: 



L + L r~ L 2 , DID " D 2 . (4) 

ft* 4 /C_4 

L + D "~ LD. (5) 

k- 5 

We assume for all reaction steps the feasibility of the reverse reaction, and that the reaction 
rates are the same for each pair of enantiomeric reactions Eqs. ([TJHD - Notice that this would 
not be the case in the presence of an external chiral polarization field, or for an internal bias 
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of non-chemical origin, such as the weak nuclear force [20|, for then k± { ^ k^. This reaction 
scheme leads to the following differential equations for the concentrations in the mean field 
limit: 

j t [L] = kx[A] + (k 2 [A] - fc_0[L] - k_ 2 [L} 2 - k- 3 [L][D] + h[A][D] - 2k 4 [L} 2 + 2k. 4 [L 2 ] 

- k 5 [L][D)+k- 5 [LD], 

j f [D] = h[A\ + (UA] - fe-i)[D] - k. 2 [D] 2 - k. 3 [D][L] + k 3 [A][L] - 2k 4 [D} 2 + 2k_ A [D 2 ] 

- k 5 [D][L]+k- 5 [LD], 

j\[A\ = -2h[A] ~ (h[A} + h[A] - fc_0([L] + [£>]) + k^ 2 ([L} 2 + [D] 2 ) + 2k^[L][D], 

j f [L 2 } = h[L] 2 - fc_ 4 [L 2 ], 
^[D 2 ] = /c 4 [-D] 2 — fc_4[D 2 ], 
^[LD] = fc 5 [L][D]-A;_5[XD]. 

The other important feature that relates our work to experimental chemical systems is the 
inclusion of both the forward ki, and the reverse reaction rates, despite the fact that 
for some rate constants one has k{ » This inclusion is necessary because these rates 
determine the thermodynamic conditions, i.e., the principle of detailed balance, microscopic 
reversibility; see for example the final paragraph of Section D 



III. OPEN FLOW SYSTEMS 



The other experimental condition that comes into play is the distinction between open 
and closed systems. Open flow is traditionally invoked for maintaining the non-equilibrium 
state. In these models, this is achieved by providing the system with a continuous input flow 
of achiral precursers A and an output flow of certain products (such as the heterodimers P 
in the classic Frank model). Openess for chemical systems is cannot be applied to recent 
experimental reports on the chemical bench (i.e., the Soai reaction). Closed on the other 
hand, refers to systems of constant mass where only energy can be exchanged with the 
exterior: the start from an initial state far from equilibrium to the final composition state is 
a real representation of a chemical reaction carried out in the laboratory. Here we consider 
reactions that take place in an open system. Open here specifically means that the system 
can exchange the substrate A with an external source, and thus we take the concentration 
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A] as a constant. This is also in keeping with the original formulation of the Frank model 



3, 111) . So we impose this condition and work out the consequences for a sequence of models 
starting from the Frank model itself. It proves convenient to use this constant concentration 
to define a simple transformation to dimensionless rates, time and concentrations. We take 
t = (k 2 [A] — k-i)t for the time parameter and \L] = ;o^z|r - [L] , etc. for the dimensionless 
concentrations. This allows us to express the above rate equations in pure dimensionless 
form: 

^-[L] = u+[L}- g[L] 2 - [L)[D] + w[D] + 2p[L 2 ] + r[LD], 
^-[D] = u+[D]-g[D] 2 -[D][L}+w[L}+2p[D 2 ]+r[LD], 

= q -[D} 2 - P[D 2 ], 
±[LD] = s[L][D] — r[LD\. 

The dimensionless parameters appearing here are: 

= k 1 [A}(k_ 3 + fcg) = fc- 2 + 2k 4 = k 3 [A] = fe-4 , , 

(k 2 [A\ - k^f ' 9 k^ + h' W (k 2 [A\ - k-i) ' P (k 2 [A]-k- 1 y U 

k-5 2/^4 k§ 
r = r , 1 ; — -, g = 77 — - , and s = — — - . (7) 

(k 2 [A] - fc_i)' (fc-S + ^s) (^-3 + ^5) 

Note we have succeeded in reducing the number of parameters from ten to seven. These latter 
seven parameters fully determine the dynamics of the kinetic scheme. We will see below 
that g plays a privileged role. Depending on whether g is greater or less than a certain 
critical value g cr u, then either the racemic or chiral solutions will be stable, respectively. 
Thus we already see the relative influence that the homo (k^) and heterodimerizations (k^) 
will have on the final outcome as well as the relative influence of the reverse autocatalytic 
(k- 2 ) and limited enantioselective steps From the denominator of g, we see that the 

chiral antagonism proceeds via two independent pathways: either through the formation of 
heterodimers (k 5 ), or via the recycling of an enantiomeric pair of monomers back to a single 
chiral monomer and the prochiral substrate (fc-3). 

To cast the equations in their final form, we next define sums and differences of the 
(dimensionless) monomer and the homodimer concentrations. Thus, we will put x = [L] + 
[D], X2 = [L 2 ] + [D 2 ], y=[L}~ [D], y 2 = [L 2 ] - [D 2 ], and P = [LD]. This then yields the 



following: 



— X = 2u + x-^(g + 1)X 2 - ^(g - l)l/ 2 + w\ + 2px2 + 2rP, 
d 

j^V = 2/(1 - w - gx) + 2py 2 , 

® Q / 2 , 2\ 

~rX2 = -Ax +y ) 



dr 
d 



-V2 



4 
<1 



PX2, 



xy-py2, 



dr" 

A P = U x 2 -y 2 )-rP. 



(8) 
(9) 
(10) 

(11) 
(12) 



dr 4' 

The fixed points of these equations correspond to the final asymptotic (as r — > oo) solutions 
of the kinetic model. Substituting x — X* + ^x(*)j y = y* + $y(t), etc. into the set of Eqs.(jH]- 
T2il where x*j2/*> denotes a fixed point solution, we obtain differential equations for the 
arbitrary perturbations 5x(t),5y(t), about the fixed point. Then, the following Jacobian 
matrix M open governs the time dependence of these perturbations to first order 0(5) in the 
fluctuations: 

j l-(g + l) X + w (l-g)y 2p 2r \ 
gy i - w - gx 2p 

p 

2 2/ 



M, 



open 



\x 



In, 1„ -p 



V 



fx 



(13) 



2 L X 

-fy -rj 

This matrix is to be evaluated on any one of the static fixed-point solutions x* i J/*> X2? 2/2 > ^* 
of the equation set Eqs.fl8l fT2|) . The associated eigenvalues indicate the stability of the specific 
static solution on which M open is evaluated. 

The eigenvalues of this Jacobian are expressed in terms of dimensionless ratios of specific 
chemical reaction rates, thus, the positivity or negativity of the individual eigenvalues can 
be ascertained straightforwardly in terms of the specific chemical parameters. This provides 
another clear connection with the experimental conditions: the relative rates of certain 
reactions will determine whether the final outcome is a chiral or racemic solution. Of course, 
here we treat the k's as variables in order to analyze all possibilities of chemical k rates. Thus 
all cases for different values and relationships between rate constants are taken into account 
in our analysis: this allows us to characterize the stable final states for reversible and quasi 
reversible, and also irreversible, synthetic chemical reactions formed by the reaction networks 
analyzed. These are the most widely studied as complex chemical systems capable of leading 



to spontaneous symmetry breaking in previous works. In practice, the rate constants can 
be determined experimentally in the laboratory, such methods are explained in detail in 
textbooks [22,|23|. 



IV. STABILITY PROPERTIES AND CRITICAL PARAMETERS: OPEN SYS- 
TEMS 

A. The Frank model 

First consider the Frank model as determined by the kinetic scheme Eqs. (fTf2f5l) . See, for 



example the definition as given in [ll|, |2l|. The concentration [A] of the achiral substrate 
is taken as constant. There is a reversible catalytic production of the monomers. The 
heterodimers [LD], formed irreversibly (&_g = 0) from the mutual inhibition step, are to be 
eliminated from the system as a kind of inactive side product. The production step Eq.flT]) 
is typically ignored (we will however consider its effect below). This situation corresponds 
then to the reaction steps Eqs. (l2f5]) . and we are to solve the equations Eqs. (IBll9l) after setting 
u = w = p = r = 0. There are four fixed points or static solutions: 

O =(x = 0,y = 0) 

R =(x = ^- g ,y = o) (14) 

O denotes the empty solution, that is, with zero chiral matter, R is the racemic solution 
with positive net total chiral matter, and Q± denote the two possible chiral solutions. From 
Eq.flSD we see that g = ^ and s = 1. 

In order to study the stability of the four possible homogeneous solutions O, R and Q±, 
we calculate the eigenvalues of the 2x2 subblock of the upper left hand corner of the 
Jacobian matrix Eq. (ll3p (after first setting w = there) evaluated at each one of the above 
four solutions Eqs. (TH|) . The eigenvalues are given by 

Ai, 2 (0) = (1,1) (15) 
Ai, 2 (i?) = (-i'T^f) ( 16 ) 
Ai, 2 (Q ± ) = (-1,=^). (17) 
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The empty solution O is always unstable. To ensure the stability in R, we need for both 
eigenvalues to be negative: Ai^-R) < 0. While the first eigenvalue \i(R) = —1 is always 
negative, we have A 2 (-R) < if and only if g > 1. Concerning the stability of the solutions Q±, 
these are always unstable if g > 1 (since \2{Q) > 0) and they become stable for < g < 1 
(since now A 2 (Q) < 0). Thus g crit = 1 [211 ] . So, the chiral solution is obtained provided that 
the rate of heterodimerization exceeds that of the reverse autocatalytic step. In the chiral 
broken solution, Q± in Eq. (j!4p . the final enantiomeric excess |ee| = \y/x\ = I =t = 1 
is always equal to unity in absolute value. These solutions are homochiral. Note that the 
irreversible heterodimer formation is mathematically equivalent to their elimination from the 
system as a continuous output, see (d) in This can be corroborated from Eqs. (TT|2][BT) . 
For elimination, then P = is constant and equal to zero in the system, and we solve the 
two equations Eqs. (|T|2|) . Instead, for irreversible formation, we set r = 0. But then Eq.(jSJ) 
for P is not independent from Eqs.flipi): we solve for (x,y) and then deduce P. The above 
steady state solutions Eq. (fT4l) and corresponding eigenvalues Eqs. (ll5HiTj) are identical for 
either situation. 

We now allow for u > and determine the effect that the direct monomer production 
step Eq.([T]) has on the solutions, and their stability properties. The four static solutions are 
now given by 



R 



= (x= 1+Vl :T 9) \ y = o) (18) 



- VA l+g IV -J 



It is easy to see that these tend to the solutions in Eq( Tl4|) in the limit as u — > 0. Here U 
denotes the unphysical solution, since the total concentration of chiral matter is negative 
X < 0. The eigenvalues are given by 



1-^1 + 4(1 + ^, 



\ li2 (R) = ( - VI +4(1 + «/)«, l + g ) ( 2 °) 



> n „ , 1 + ^l + ±g(-l + g + ±g 2 u) -1 + y/1 + Ag(-1 + g + 4g*u) ^ 
Ai, 2 (Q±) - ( Tg , Yg )■ (21) 



In the u —>■ limit these eigenvalues tend those listed in Eqs. fll5HiTj) . Since Ai j2 (f/) > 0, the 
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unphysical solution is always unstable. Note also that Xi(R) < and Xi(Q) < are always 
negative whereas \2(R) > and X^Q) < for g < g crit where g cr i t = g^(Vl + 16m — 1). 
Note that g cr it{u) < 1 for all u > 0. For small u we can write g CT u = 1 — 4u; while for large 
u, g cr it —> 2m 1/ 2 ■ Thus the monomer production step tends to racemize the system, lowers 
the final ee to values strictly less than unity: 



ee = ±^/l-i^, (22) 

and drives the value of g cr u below one. The monomer production step thus reduces the 
range of g for which stable mirror symmetry breaking can occur, and the chiral solutions 
are no longer homo chiral. 

B. Frank with reversible heterodimer formation 

What happens if the heterodimers are not removed from the system, but are allowed 
to remain evolving dynamically in concert with the monomers? This situation corresponds 
to all three reaction steps Eqs. (nni2l[5|) . and we solve the three equations Eqs. (l8f9][T2"]) after 
setting w = p = q = and s — 1. To keep the algebra manageable, we also will set u = 0. 
There are four static solutions: 

O =(P = 0, X = 0,y = 0) 

R ^(P=^, x = ^ = 0) (23) 

Note that the final heterodimer concentration P is zero in the chiral states Q±. 

In order to study the stability of the four possible homogeneous solutions O, R and Q±, 
we calculate the eigenvalues of the 3x3 array obtained from Eq. tflBl after deleting the 3rd 
and 4rth rows and columns. The eigenvalues corresponding to these solutions are given by 

A 1)2)3 (0) = (1,1, -r) (24) 
2 + g(l + r) + y/4 + g 2 (-l + r) 2 + 4g(l + r) 

2g 



Al,2,3(-R) 



Xi^iQz 



-2 - g(l + r) + yj + g 2 (-l + r) 2 + Agjl + r) 

2g } 

/_ l+g(-l + r) + y/l + 2g(-l + r) + ^ 2 (1 + r) 2 
V ' 2g 



(25) 
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-1 + g(l - r) + + g(-l + r)) 2 + Ag 2 r 
2~9 



(26) 



As Ai ) 2(0) > 0, the empty state is always unstable. An inequality analysis shows that both 
\2(R) < and X-s^R) < for all r > and g > 0. Since Xi(R) = —1 this demonstrates that 
the racemic state R is always stable. As an independent check, we also verify that X$(Q) > 
is positive for all r > and g > 0, so the chiral solutions Q± are always unstable. There 
is no stable mirror symmetry broken solution when the heterodimers (formed through the 
crucial mutual inhibition step) are included reversibly in the system, see Table [H 



C. Frank with reversible homo- and heterodimer formation 

If heterodimers can form, it is certainly reasonable to expect the same for homodimers. 
This corresponds to the reaction steps Eqs. ([Tll2Tfl5l) . and we must solve all five equations 
Eqs.( l8|9|10lllllll2l) after setting w = and s — 1. To help maintain algebraic control, we 
also will set u = 0. There are four static solutions: 

O = (P = 0, X2 = 0,y 2 = o, x = o,y = 0) 

It is important to point out that g = fc ~ 2 ^ 2fc4 = ^ + q, so that g > q. Note that the 
final heterodimer concentration P is zero in the chiral states Q±, but there a net positive 
concentration \i > for the homodimers. 

In order to study the stability of the four possible homogeneous solutions O, R and Q±, 
we proceed to calculate the eigenvalues of the full 5x5 array in Eq. (fr3l) . As it turns out, we 
are unable to obtain analytic closed-form expressions for the eigenvalues of the racemic state 
R, but we can do so for both the empty O and the chiral solutions Q±, which is sufficient 
for our purposes. These are given by: 

Ai,2,3,4, 5 (0) = (1,1, -p,-p,-r) (28) 

9 + gp + q-pq+ yj-4p{g -q) 2 + (g + gp + q- pq) 2 



Al,2,3,4,5(Q±) = i~Pi 



2(9 - q) 



-g(l+p) + (-l + p)g+ y/-4p(g - g) 2 + (g + gp + q-pq) 2 

2(9 ~ q) 
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1 + q + g(-l + r) - qr + a/4(# - q) 2 r + (1 + q + g(-l + r) - gr) 2 

% - g) 



-1 + 5- - g - #r + qr + y/4(g - q) 2 r + (1 + q + #(-1 + r) - gr) 2 \ 

> (29) 

As before, the empty solution is always unstable. As for the chiral solution, an inequality 
analysis shows that \\ y 2,3,i{Q) < are all negative, whereas X^(Q) > 0. Hence the chiral 
asymmetric solutions are unstable, and thus by logical deduction, the racemic state R must 
therefore always be stable. We point out however, that if the heterodimers are formed 
irreversibly (r = 0), then a stable homochiral outcome is again possible whenever g < 1. 

To summarize up to this point, the sequence of Frank type models analyzed above indi- 
cates that the elimination of the heterodimers from the system (when [A] is constant), or 
their irreversible formation, is crucial in order that symmetry breaking be possible and that 
the chirally asymmetric states be stable. Allowing for homodimer formation does not alter 
the stability properties. See also the first four rows of Table HI 

D. Limited enantioselectivity 

Whereas the forward reactions in Eq.([2]) represent the autocatalytic capabilities of each 
enantiomer, the forward reaction in Eq. ([3]) accounts for the fact that an enantiomer can also 
catalyze the production of its chiral partner. This scenario, corresponding to the reaction 



steps Eqs. (TT|2][5]) . is termed limited enantioselectivity [13j, |24|. As there is no dimerization 
in this model, we are to solve Eqs.( l8ll9l) in which we set \i = P = 2/2 = 0. In this case, the 
four static solutions are (note we also set u = to help simplify the algebra): 

O =( x = 0,y = 0) 

R =(x = ^,y = o) (30) 

O denotes the empty solution, that is, with zero chiral matter, R is the racemic solution with 
positive net total chiral matter, and Q± denote the two possible chiral solutions. Because 
there is no dimerization in this model, we can set = k§ = in Eq.([6]) and so we find 

In order to study the stability of the four possible homogeneous solutions O, R and Q±, we 
calculate the eigenvalues of the 2x2 subblock of the upper left hand corner of the Jacobian 
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matrix Eq. (fT3"|) evaluated at each one of these four possible solutions. The eigenvalues are 
given by 

A lj2 (0) = (l-w,l + w) (31) 
\ (b\ t 1 1-w - g - 3wg 

Ai, 2 (i2) = {-1-w, ^— ) ( 32 ) 



Ai, 2 (g±) 



-l+w + 2wg- - w) 2 + g 2 (A + 8w - 8w 2 ) - Ag(l - w) 

2^ 

-l+w + 2wg + - w) 2 + g 2 (A + 8w - 8w 2 ) - Ag(l - w) 

2~g 



(33) 



The empty solution O is always unstable. To ensure the stability in R, we need for both 
eigenvalues to be negative: Xi^iR) < 0. While the first eigenvalue is always negative, we 
have A 2 (i?) < if and only if g > jt^- Concerning the stability of the solutions Q±, these 
are always unstable if g > 1 and they become stable for < g < < 1. When this holds, 
the final enantiomeric excess ee will be given by 




ee = ± h_ (34) 

and so < |ee| < 1, which means that the chiral solution is not homochiral. This symmetry 
breaking can be represented via a corresponding bifurcation diagram, see Fig [U When 
g > g C rit only the racemic state is stable, and ee = 0, but when g < g crit two stable chirally 
asymmetric states are equally likely, which one the system actually chooses is a random 
event. The upper and lower branches of the bifurcation are given by plotting the values of 
ee, Eq.pl]) function of g holding w fixed. 

In summary, there is a critical value of the parameter, 

1 + 6W 

that uniquely determines the outcome of the reaction scheme. For g > 1^™ W the only stable 
solution is R. In this case y = i.e. [L] = [D], i.e. we have a racemic solution. Note 
that if the decay k-± of the enantiomers into achiral matter can be neglected with respect 
to the rate of autocatalytic amplification, then w is well approximated by the ratio ^3/^2, 
independent of the concentration of achiral matter [A], and so the critical value g cr u and 
final solution are controlled by the competition between the rates of the forward reactions 
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Eqs.(2.2) and (2.3). That is, the relative rate of autocatalysis versus limited enantioselec- 
tivity, in which the catalytic effect of each enantiomer leads to the formation of both L 
and D products. On thermodynamic grounds however, the symmetry breaking condition 



g < g cr it, determined from the stability analysis, cannot be achieved [241 ] . Enantiomers are 
thermodynamically identical species and therefore, at the final state must fulfill the condition 
= = which is a consequence of the principle of detailed balance 25|]. Therefore, 
the necessary condition g < 1 (i.e., k_ 3 > k_ 2 ) is incompatible with k 2 > k 3 , and g will 
always be greater than the critical value. The question then arises if it is possible or not to 
find more complex non-linear reaction networks where the mathematical condition for sym- 
metry breaking, as determined by the stability analysis, is compatible with thermodynamic 
chemical constraints. 



bifurcation diagram 



a-b 




unstable 
racemic 
solution 



racemic 
solution 



8568 Of 



8569 85695 0.B57 85705 8571 85715 0.8572 
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FIG. 1: Bifurcation diagram, limited enantioselectivity. For g > pfg^ , the unique solution is 
y = 0, i.e. [L] = [D], which corresponds to a racemic solution. Reducing g below the critical value, 
the system undergoes a pitchfork bifurcation, leading to chiral solutions, i.e. [L] ^ [D]. This is 
illustrated here for g cr it = 0.85714 (w = 0.04). These chiral solutions however are not homochiral, 
as they yield an enantiomeric excess ee strictly less than unity, provided that ^ 0; see Eq. fpP 



E. Limited enantioselectivity plus dimers 

It is natural to ask if including dimer formation can alter the symmetry breaking condition 
9 < gcrit and get around the thermodynamic constraint. To this end, and for more chemical 
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TABLE I: Selected reaction schemes treated as open systems, and the corresponding condition 
9 < Qcrit that must be satisfied for obtaining a stable chiral solution. 



Reaction model (open systems: [A] = constant) g 



Qcrit 



9 < 9crit 



Frank (with fa 



0) a 



fc-2 
&5 



k-2 < fa 



Frank (with fa > 0, k-i > 0) 



Frank + reversible heterodimers 



fc-2 
k 5 



Frank + reversible homo- and heterodimers e 



fc-2+2fc 4 
&5 



2A: 4 
k 5 



Limited enantioselectivity 9 



k-2 
k- 3 



l—w h 



l+3w 



fc-2 , l—w 

ks l+3w 



Lim. enantiosel. plus heterodimers 



fc_3+fc5 



1— U! I 



l+3w 



fc-2 ^ 1-tu 
fc-3 1+3TO 



Lim. enantiosel. plus homodimers 



fc_2+2fc 4 
fc-3 



l+q— w+Zqw j 
l+3w 



fc-2 ^ 1-w 
fc-3 l+3u> 



Lim. enantiosel. plus homo- and heterodimers k k 1+9 s ±+3^, 9W+sw k 



fc-2 , 1-W 
fc-3 l+3w 



"Reaction steps Eqs. (I2I5P ; see Sec IIV Al The system is fed by an input of the achiral substrate A, and the 



11 



2l|. 



output consists of the LD- heterodimers formed in the mutual inhibition step Eq.([5]); see e.g., Ref. 2, 
^Reaction steps Eqs.(HGEl); see Sec lTVAl Here, u = r^^Z^ 

c Reaction steps Eqs. (I2I5P ; see Sec IIV Bl The heterodimer concentration is time dependent, in contrast to 

the original Frank model. 
d 9crit = 1 — s, but in this model k 3 — k-3 — 0, so that s — fa/k§ = 1. 

e Reaction steps Eqs. (|2l4|5|) : see Sec IIV Cl Both homo and heterodimers are formed, and are not removed 
from the system. 

^In this model, g > q = 2fc 4 /fc5, so no symmetry breaking is possible. 

^Reaction steps Eas. pTg]) , see Sec HVD] 

''Recall w = fa/k2 for large constant concentrations [A] . 

'Reaction steps Eqs. (I2I3I5|) . see Sec IIV El In this model s = fc 5 /(fc„ 3 + fa) < 1. 

^Reaction steps Eqs. (|2l3l4p , see Sec IIV El In this model, since fcg = 0, then q = 2fc 4 /fc_3. 

^Reaction steps Eqs. (|2"llo]) . see Sec IIV El Here q = 2k A /(k- 3 + k 5 ) and s = fa/(k- 3 + fc 5 ) < 1 for fc_ 3 > 0. 
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realism, we include the dimerizations. We begin with just the hetero-dimerization step, and 
the static solutions work out to be given by 



O = (P = 0,y = 0, X = 0) (36) 
s(l +wf _ _ 2(1 + w) 
r(l + g — s) z 1 + g — s 



^ = ( p= , : ^2 »y = o.x= , , i J (37) 



The eigenvalues for the empty and racemic solution can be found in closed form. They 
are given by 

\i, 2 , 3 {0) = {-r,l-w,l + w). (39) 



and 



/ , 2g(l + w) 



1 + g-s 1 
1 + r + s — rs + w + sw + g(l + r + w) 
2{l + g-s) 



y — 4r(l + g — s) 2 (l + w) + (1 + r + s — rs + w + sw + g(l + r + w)) 2 
+ 2(l + g-s) ' 

r(-l - g + s) - (1 + g + s)(l + w) 

2(l + g-3) 

a/— 4r(l + (7 — s) 2 (l + w) + (1 + r + s — rs + w + sw + g(l + r + w)) 2 \ 
+ 2(l+g-s) J' ( 0) 

Note that here, g = k h ~+ k5 ■ The empty solution O is always unstable. As for the racemic 
solution R } we note that \\(R) < provided that g > • The remaining two 

eigenvalues A2,3(-R) < for r > 0, w > 0, 1 > s > and for all g > 0. So the racemic state 
is stable when g > (7 crit = and unstable when ^ < (7 cr jj. Clearly, when the racemic 

state is unstable, we infer that chiral state is the only stable solution. 

Note that the expression for the critical g in this model is just the g cr u for the limited 
enantioselectivity model treated above times the additional factor (1 — s), see Table [H In 
other words, the inclusion of heterodimer formation and dissociation seems to reduce the 
value of g CT iti appearing to making chiral symmetry breaking relatively more difficult to 
achieve. But since s = k ^ k5 , the crucial inequality g < g crit actually reduces algebraically 
and simplifies to yield the one obtained for the limited enantioselectivity model without 
heterodimers. 
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For completeness, we also have worked out the fixed point solutions and their stability 
properties first, when only the homodimers, and then when both the homo-and heterodimers 
are included. The results of these latter two cases are summarized in the seventh and eighth 
rows of Table [U It is interesting to point out that although the expressions for the critical 
g's depend on the homo and heterodimer rates explicitly and in a non-trivial way (see the 
third column), the final inequalities g < g cr u determining the stable chiral outcomes are 
the same as for limited enantioselectivity alone, irrespective of whether the homo and/or 
heterodimers are included. Compare the fifth and seventh rows of the Table, as well as the 
sixth and eighth rows. The dependence on k± or q simply cancels out. The inclusion of 
dimer formation does not affect the symmetry breaking in this model. Since the condition 
■p- = i~~ = -r 2 - continues to hold, only the racemic state is stable. 



V. CLOSED REACTION SYSTEMS 



For closed systems there is no flow of material into or out of the system. Since [A] is 
not constant in this situation, we cannot use it to rescale the time or the concentrations. 
We can however take r = kit for the time and [L] = 5 [L] , etc. for the dimensionless 
concentrations. This allows us to express the rate equations in Sec [TT] in the following 
dimensionless form: 

±[L] = [A] - u[L] + h[A][L) - g[Lf - [L][D] + r[A][D) + 2u[L 2 ] + p[LD], 
j-[D] = [A] - u[D) + h[A)[D] - g[Df - [D][L] + r[A][L] + 2u[D 2 ] + p[LD), 

= q -[Lf-uj[L 2 l 
^[D2} = q ~[D] 2 - u>[D 2 ], 
±[LD] = s[L\[D\-p[LD\. 

These are subject to the constraint [A] = C — [L] — [D] — 2[L 2 ] — 2[D 2 ] — 2[LD], where C is 
a constant. The parameters appearing here are 

k-i k^ 2 + 2/^4 k 2 ks k-5 A;_4 2k^ 

u = i— > 9 = ~i n—i h = i n-i r = i —r,P = -^,u = -^i(l 



h ' A;_ 3 + k 5 ' A;_ 3 + k 5 ' /c_ 3 + k 5 ' kx k\ ' A;_ 3 + k 



5 



and s = — — 5 , . . (41) 
(k-s + h) 
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For closed systems, the number of independent parameters reduces from ten to eight. 

As before, we find it convenient to define the sums and differences of concentrations: 
x =[L] + [D},y= [L] - [D], X 2 = [L 2 ] + [D 2 ] y 2 = [L 2 ] - [D 2 ], and put P = [LD]. This 
then yields the following: 



= 2iA] + ((h + r)[A]-u)x-\(g + l)x 2 -l(g-l)y 2 + ^X2 + 2pP, (42) 



j^V = y{{h - r)[A) - gx - u) 



2w y 2 , 



d 

dr 
d 



X2 = |(x 2 + y 2 ) -wx2, 
d -P = -Ax 2 - y 2 ) - pp. 



(43) 
(44) 
(45) 
(46) 



dr 4' 

In these variables, the constant mass constraint reads [A] = C — x ~ %X2 — 2P. 

For the complete closed model, that is, the reaction steps in Eqs. (JTJ{5j) , the 5x5 Jacobian 
matrix M c i ose d for the linearized fluctuation equations is given by 

/ thj-*-\(n-i^ r .-iP\-0-i.„\ {l-g)y 2w-A-2(h+r)x 2p-A~2(h+r)x « 



M, 



closed 



(h+r)(C-2 X2 -2P)~(2+u) 
-(g+l+2(h+r)) X 



\ 



+h—r)y 


{h~r)(C-2 X 2-2P) 


-2(h-r)y 


2u> 


-2(h-r 




-(g+h~r)x~u 








lx 




—UJ 








\v 


fx 





— OJ 





fx 










-p 



(47) 



/ 



This incorporates the constraint equation directly. This matrix is to be evaluated on any 
of the static fixed point solutions x*i 2/*; X25 V2 an d P* °f the equation set Eqs. (l42ti46T) . The 
constant C is related to the total initial concentration Q as follows: C = k ~ 3 k ^ k5 Q. The 
eigenvalues of M c i ose d indicate the stability of the static solution on which this matrix is 
evaluated. Unlike the open flow case however, it is much more difficult to obtain manageable 
analytic expressions for these eigenvalues, and the variety of reaction schemes we can treat 
analytically in this manner is limited. 
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A. Closed Frank with irreversible heterodimers: fc_5 = 



Recall that the Frank model in an open flow system and with irreversible formation of 
heterodimers does lead to stable steady-state chiral solutions. We now enclose it in a box, 
so that [A] is no longer constant. In this situation we are to solve the Eqs.( l42|l4lIf4"6T) where 
from the parameter list in Eq (j4"T]) we have u > 0, g = k^/ks, h = A^/^s, r = u = q = 0, 
s — 1 and for fc_ 5 = 0, we must set p = 0. The mass constraint reads [A] = C — \ ~~ 2P. 
There is a static racemic R and two static chiral Q± solutions: 

R =(P = q,y = 0, X = 0) 
< ?±s (i»=2g=, s = ±;, x =-j). 

Note that the chiral solutions are unphysical since they imply % < when u > 0. We 
can attempt to get some information regarding stability by evaluating the 3x3 array ob- 
tained from Eq. (j4"Tl) after deleting the 3rd and 4rth rows and columns. The eigenvalues 
corresponding to the racemic solution are 

A 1)2>3 (i2) = (0,-2-u,-w), (49) 

indicating that R is marginally stable (marginal because of the zero entry). We are unable 
to calculate the eigenvalues for the unphysical chiral states in simple analytic form. For 
u = 0, the only steady state is the racemic one R, (there is no steady chiral solution) and 
its eigenvalue is (0,-2,0). 

n 

This is important, because Model 1 of Rivera Islas et al [8J, is actually a special case of 
this, in which u — g — 0, and we consider this next. 



B. Rivera Islas et. al. Model 1 



This model 



corresponds to our reaction steps Eqs (fT|2|5|) and after setting k_ x = k_ 2 = 
fc_5 = 0. Note this model implies irreversible reactions. This is a closed system, we we solve 
again the Eqs. fl42f43ir4^|) where from the parameter list in Eq fl4ip . we now have u = 0, g = 0, 
h = k 2 /k 5 , r = u = q = 0, s — 1 and for fc_ 5 = 0, we must set p = 0. The constraint is 
[A] = C ' — x~ 2-P '. We next look for the static solutions of the set of equations Eqs. (l42f43f46l) . 
The racemic solution y = implies that x = an d P = C/2. Thus, all the net matter 
ends up as heterodimers. On the other hand, the chiral solution has y ^ 0. For this case, 
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the static solution of Eqs. (j42j43|4U|) yields x 2 = U 2 an d C — IP = X - The physical solution 
corresponds to a positive x > 0, so that y = ±x an d P = ^-rp > 0. In this case, the total 
matter is distributed among monomers and dimers. Note however we cannot independently 
solve for the final heterodimer and monomer concentration in the chiral phase. The static 
solutions are summarized as follows: 

R =(P = q,y = lX = 0) 

Q ± = (p=^2C,y = ± X ). 

Using M dosed to compute the eigenvalues for this case, we find that 

K^(R)isias = (-2,0,0). (51) 

Since \2,3{R) T , = 0, we cannot say conclusively if this state is stable or unstable. This 
is in accord with the chemical fact that due to its irreversibility, this reaction scheme leads 
to a kinetic controlled outcome of the reaction products. Evaluating Eq. (1471) on this chiral 
solution and calculating the eigenvalues yields 

Xi,2, 3 (Q)isi as = ( Q ,-X,-2-h X ). (52) 

This is as much as we can say regarding linear stability analysis; there are no free adjustable 
parameters that can induce a change in sign in any of the eigenvalues. The simple stability 
analysis is inconclusive. Nevertheless, numerical simulations carried out in [8( indicate stable 
chirally asymmetric solutions for a range of h. Thus, a higher order stability analysis might 
be called for. 



C. Limited enantioselectivity revisited: the closed system 

We can obtain the static solutions and their associated eigenvalues exactly for limited 
enantioselectivity in a closed system, the model introduced in Subsection IIVDI This is 
useful because it provides an exact point of comparison between a specific set of reactions 
in both open and closed systems. So we consider the reaction set Eqs.( fTH2ll3l) . but now in 
a closed system. We must solve Eqs. fl42f43l) in which we set X2 — P — V2 — 0- The mass 
constraint in this case reads [A] — C — x- There are four static solutions (note we also set 
u = to simplify somewhat the algebra): 

tj= (v = 2-C(h + r) + y/4C(l + g + 2(h + r)) + (2 - C(h + r)f 

{V ' X l + g + 2{h + r) 
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-2 + C(h + r) + V4C(l+ g + 2(Hr)) + (2-gCT 
l2/ ' X l + # + 2(/i + r) 

n , , l *Cg(g + h-r) + C 2 (h - r)((-l + ^ + r + Zg~r) C(h-r) 

Q± = iy = ± V (g-l)(g + h-ry ' X = JTh^r^ 

U denotes the unphysical solution, as this implies a negative total chiral matter, R is the 
racemic solution with positive net total chiral matter, and Q± denote the two possible chiral 
solutions. Note from Eq. (j4ip we have g = h = r = 

The stability of the four possible homogeneous solutions U, R and Q±, is determined 
from considering the eigenvalues of the 2x2 subblock of the upper left hand corner of the 
Jacobian matrix M dosed Eq. (j4"7"|) (and after setting \2 = P = as well as u = 0) evaluated 
at each one of these four possible solutions. The eigenvalues are given by 

Ai, 2 (C0 = (C(h - r) - ^Jg±^L (_ 2 + C(h + r) - ^4 + C\h + rf + 4C(1 + <? + ft + r)) , 
^4 + C 2 (/i + r) 2 + 4C(1 + g + /i + r)) (56) 



Ai, 2 (i2) = (C(h-r)- - r^Ar^k (-2 + C(h + r) + ^A + C^h + ry + AC(l + g + h + 
V 1 + a + In + 2r V 



1 + + 2/i + 2r 



r 



- v '4 + C 2 (/i + r) 2 + 4C(l + + /i + r)J (57) 

and 

Ai, 2 (Q±) = ^T—^ Mg(-l + Cr) + (-h + r)(2 + C(l + h + r)) 

2{g + n — r) 



+ y/((g(2 - 2Cr) + {h - r){2 + C{1 + h + r))f 
r)(4g 2 - C(h - r) 2 + g(h - r)(4 
(2 5 (-l + Cr) + (-/i + r)(2 + C(l + /i + r)) 



+ 4C(# + /i - r)(4# 2 - C(/i - r) 2 + - r)(4 + C(/i + 3r))))), 
1 



2(<7 + /i - r) 

- y((</(2 - 2Cr) + (h-r)(2 + C(l + h + r))) 2 

+ 4C(# + h- r){4g 2 - C(h - rf + g{h - r)(4 + C(h + 3r))))) 

As A 2 (£7) > 0, the unphysical solution U is therefore always unstable. To ensure the stability 
in R, we need for both eigenvalues to be negative: Ai )2 (i2) < 0. While the second eigenvalue 
is always negative, we have Xi(R) < if and only if g > g cr u, where 

gent = -g(/i - r) (A + C(h + 3r) - ^16 + C 2 (h + 3r) 2 + 8C(2 + h + 3r)) . (58) 
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A careful expansion of this shows that in the limit of large C, then the critical value of g 

_ 1-r/h _ 1 - h/k 2 _ 1 - w _ open 

- x + 3r/h - x + 3h/h ~ 1 + 3w ~ 9cnt , ^y; 

asymptotically approaches the critical g for the same set of reactions in an open system, 
Eq.fl35D. 

Concerning the stability of the solutions Q±, these are always unstable if g > 1 and they 
become stable Xi,2(Q±) < for < g < g cr i t < 1. However, the thermodynamic condition 
mentioned above still holds and implies that the racemic state is the only stable outcome. 



VI. CONCLUDING REMARKS 



We have focussed attention on the steady state solutions and their dynamic stability 



properties in variants of both the Frank anc 

r 

symmetry breaking. In the Frank model [7, 



imited enantioselectivity models of mirror 



the mutual inhibition occurs through 
the formation of heterodimers composed of the two chiral monomers, whereas in limited 



enantioselectivity 13|, |24j , the needed chiral antagonism occurs though a monomer recycling 
reaction. We have unified both models into one encompassing reaction scheme and have 
studied the steady state solutions and their stability properties in a sequence of reaction 
schemes that start from the paradigmatic Frank model to limited enantioselectivity, in which 
we study the consequences of including the formation of (reversible) hetero- and homo-dimers 
as dynamic variable concentrations in both. The general conclusion we can draw from this is 
that the inclusion of the variable heterodimer concentrations leads to a final stable racemic 
state in the sequence of Frank-type models defined above. This is important because we 
recall that in Frank's original formulation, the heterodimers are supposed to be removed 
from the system, i.e., their formation is irreversible. Including the homodimers does not 
affect this conclusion. 

In the case of the limited enantioselectivity model, no dimer formation is originally con- 
templated, and mirror symmetry breaking is mathematically possible because the needed 
inhibition is provided by the reverse symmetric catalysis step [24 |. In spite of this, the 
symmetry breaking condition g < g cr u, determined from the pure stability analysis, cannot 



be achieved on thermodynamic grounds [24|, as has also been pointed out recently 26]. 



Including the homo- and/or the heterodimer formation to limited enantioselectivity does 
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not alter the crucial inequality, nor does the presence of the dimers change the thermody- 
namic constraint. Therefore an open question is if a different chemical scenario leading to 
collective phenomena exists, such as a second order phase transition, that can give rise to 
the condition g < g cr u- In the context of spatially extended polymerization systems, it has 



been argued before that one s 
order phase transitions 



nou 



27 



28 



29]. 



d expect domain formation similar to the case of second 



The reaction networks treated here are limited to reactions involving only monomers and 
dimers, whereas biological chirality of living systems involves large macromolecules that are 
probably the result of polimerization reactions. In this vein, Sandars recently introduced 
a model in which the detailed polymerization process and enantiomeric cross-inhibition 
are taken into account, its basic features are explored numerically, but without including 
spatial extent, chiral bias or noise jsO]. Brandenburg and coworkers have analyzed further 
properties of Sandars' model and have proposed a truncated version including chiral bias [3l| , 



and have studied this reduction with spatial extent and coupling to a turbulent advection 
velocity [27]. Gleiser and Thorarinson analyze the reduced Sandars' model with spatial 
extent and coupling to an external white noise [32| and in [28( , Gleiser considers the reduced 
chiral biased model with external noise. In addition to Sandars, both Wattis and Coveney 



331 ] and Saito and Hyuga 34| have introduced polymerization models that can give rise to 



homochiral states. The latter one differs from Sandars' in allowing for reversibility in all the 
steps. As reported in the review article (d) in 9], all these polymerization models derive 
from Frank's model by adding polymerization reactions. Thus, in spite of the simplicity of 
the original Frank model ignoring as it does the polymerization process, it continues to 
serve as a type of "Ising model" for chiral symmetry breaking. 
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